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We introduce a method for accurate quantum chemical calculations based on a simple variational 
wave function, defined by a single geminal that couples all the electrons into singlet pairs, combined 
with a real space correlation factor. The method uses a constrained variational optimization, based 
on an expansion of the geminal in terms of molecular orbitals. It is shown that the most relevant non- 
dynamical correlations are correctly reproduced once an appropriate number n of molecular orbitals 
is considered. The value of n is determined by requiring that, in the atomization limit, the atoms are 
described by Hartree-Fock Slater determinants with Jastrow correlations. The energetics, as well 
as other physical and chemical properties, are then given by an efficient variational approach based 
on standard quantum Monte Carlo techniques. We test this method on a set of homonuclear (Be 2 , 
B 2 , C 2 , N 2 , 2 , and F 2 ) and heteronuclear (LiF, and CN) dimers for which strong non-dynamical 
correlations and/or weak van der Waals interactions are present. 



I. INTRODUCTION 

Already in the early stages of quantum mechanics, L. 
Pauling introduced the so called resonating valence bond 
(RVB) theory of the chemical bond[l[, starting from the 
simple consideration that a spin singlet can be formed be- 
tween any two valence electrons belonging to neighboring 
atoms. In this scheme, the ground state wave function 
of a molecule, such as benzene, can lower the energy by 
allowing the resonance among all possible valence bond 
configurations that can be drawn by linking the positions 
of two atoms (e.g. the Kekule and Dewar configurations 
in the benzene molecule). However, its application was 
limited, since the number of bonds were growing expo- 
nentially with the number of atoms. As a consequence, 
the powerful language of molecular orbitals (MO's) ap- 
plied to Hartree-Fock (HF) and post HF methods be- 
came popular. Nonetheless, quite recently, the interest in 
RVB wave functions has been strongly revived. Indeed, 
soon after the discovery of the High-Tc superconductors, 
P. W. Anderson realized that a single determinant wave 
function combined with a suitable real space correlation 
term - henceforth referred to as 'the Jastrow factor' - 
could be used to represent a complex RVB state. [2[ In 
this new ansatz a crucial ingredient is the form of the 
determinantal part of the wave function, that is required 
to be a singlet state with total spin 5 = 0. This picture, 
aimed at explaining the High-Tc superconductivity, rep- 
resents also a very efficient numerical implementation of 
the original RVB idea^soon reconsidered in this form for 
lattice models 0, 0, S E , 

and then in realistic simula- 
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tions of atoms and small molecules 

BStMtlil- Though 
the Anderson's RVB wave function has been originally 
defined just for singlet states, the same concept can be 
applied to electronic systems with arbitrary spin S > 0, 
with the inclusion of unpaired orbitals. This is a very 
important generalization in order to describe polarized 
compounds, like the transition element compounds which 
show high-spin configurations in their low-lying energy 
states. In the actual RVB description of realistic sys- 
tems, it is necessary to resort to standard quantum Monte 
Carlo (QMC) methods [t3| in order to compute the vari- 
ational expectation values of the energy and correlation 
functions 0, H, [H . 

In this article, we propose an extension of the RVB pic- 
ture that is based on a MO expansion of the singlet va- 
lence bond pairs defining the wave function. This ansatz 
yields a correlation consistent RVB representation by 
means of a constrained energy minimization which keeps 
the number of MO's fixed while stretching the bond. By 
setting this number to a value such that a Jastrow cor- 
related HF wave function is recovered in the atomization 
limit, we obtain, with a single determinant, a remark- 
ably accurate description of the bond, even when strong 
non-dynamical correlations are present in the system. In 
this paper we illustrate the method and test it on a set 
of dimers composed by first row atoms and on selected 
small molecules belonging to the so-called 'Gl set' (sec 
Ref. [Dj])) often used to test new theoretical methods. 
The approach described in this work has also been ap- 
plied to the study of the controversial ground state of the 
iron dimer[l5l|. 

In the following, we describe the RVB wave function 
and our extension, and we show test results on various 
homonuclear and heteronuclear dimers (Be 2 , B 2 , C 2 , N 2 , 
2 , F 2 , LiF, and CN). In Appendix |A1 we describe the 
constrained minimization of the molecular orbital expan- 
sion of the trial wave function. In Appendix IB"! we present 
a systematic study of the variational energies obtained 
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with the RVB wave function, as a function of the num- 
ber of molecular orbitals (Appendix IB 1[) and the size of 
the atomic basis set used (Appendix IB 2|l . 



II. VARIATIONAL METHOD 
A. General description of the wave function 

The fundamental ingredient of our variational method 
is an N-electron RVB wave function, called JAGP since it 
is the product of a Jastrow factor J, and a determinantal 
part which is an antisymmetrized geminal power (AGP), 
previously introduced in Refs. [8, 9] (^jagp = J&n)- 
Below we shall describe this wave function. 

In the case of N electrons with N-\ up spins (N± down 
spins), where for simplicity we take N-\ > iVj., we can de- 
scribe a pure spin state with total spin S = \Nf — iVj_j/2 
and maximum spin projection S 1 ' ' = S by means of the 
antisymmetrized product of N± singlet pairs and 2S un- 
paired orbitals corresponding to the remaining spin-up 
electrons. Hence the determinantal part reads 

.a n n w 

i=l j=N l + l 

with A the antisymmetrization operator, R = 

|r}, • • • , r]y T , r{, • • • , rj^ | the 37V-dimensional vector of 

coordinates, (f>(r,f') — tpif ,r) a symmetric orbital func- 
tion describing the singlet pairs, and <f>j(r) the unpaired 
orbitals. It can be shown that the wave function in 
Eq. ([T]) can be rewritten in terms of a single determinant 
(see Ref. [§| and references therein). 

$n(R) has definite total spin. We also impose all pos- 
sible symmetries to be satisfied by Eq. (fT]), such as an- 
gular momentum and spatial reflections. 

Similar constructions with definite spin can be done, 
by allowing also triplet pairing between the 2S unpaired 
electrons. Since this involves a bit more complicated al- 
gebra like the use of pfaffiajis [10. Il l| . we do not consider 
it here. 

The Jastrow factor takes into account the electronic 
correlation between two electrons and is conventionally 
split into a homogeneous interaction Ji depending on the 
relative distance between two electrons (i.e. a two-body 
term), and a non homogeneous contribution depending 
also on the electron-ion distance, included in the one- 
body J\, three-body J3 and four-body J4 terms. J\ is 
a single particle function which is important to compen- 
sate the change in the one particle density induced by 
J2, J3 and J4, as well as to satisfy the electron- ion cusp 
conditions. The one- and two-body terms J± and J2 are 
defined by the following equations: 

.h = exp -{2Z a f^u{Z\l\ ia ) + J29?X J al (r t )] , (2) 

ia ial 



and 

J 2 =exp[^u(r ii )], (3) 

i<j 

where i,j are indexes running over the electrons, and I 
runs over different single particle orbitals xti centered 
on the atomic center a. Ti a and ry denote electron- 
ion and electron-electron distances respectively. The 
corresponding cusp conditions are fixed by the function 
u(r) = F[l - exp(-r/F)]/2 (see e.g. Ref. pj]). at and 
F are optimizable variational parameters. 

The three- and four-body Jastrow J3J4 are given by: 

J 3 J 4 (R) = exp \T f^' rj) j , (4) 

with f(r,f), being a two-electron coordinate function 
that can be expanded into the same single-particle basis 
used for J\\ 

m,?i) = E at xti(n)xLm (5) 

ablm 

with g?^ optimizable parameters. Three-body (electron 
ion electron) correlations are described by the diago- 
nal matrix elements g aa , whereas four-body correlations 
(electron ion electron ion) are described by matrix ele- 
ments with a =/= b. 

The exhaustive and complete expression of the Jastrow 
factor J(R) = J 1 (R)J 2 (R)J 3 (R)J4(R) that we adopt 
in this work allows to take into account not only weak 
electron-electron interactions of the van der Waals (vdW) 
type, but it is also extremely effective for suppressing 
higher energy configurations with overlapping valence 
bonds, which otherwise lead to a too large electron den- 
sity around an atom. 

As any functions of two coordinates, also the pairing 
function <j) in Eq. ([T]) can be expanded in terms of single 
particle orbitals. We can thus write 

n-25 

Hry)= J2 WiWM?), (6) 

3=1 

where n is large enough, and {<fij} is an orthogonal single 
particle basis set, which reaches its complete basis set 
limit (CBS) for n — > 00. Notice that, in these notations, 
we assume that the 2S unpaired orbitals <f>j of Eq. (Q]) 
correspond to the indexes n — 2S + 1 < j < n in Eq. ©. 

The single particle orbitals 4>j can be conveniently cho- 
sen as the MO's obtained with a conventional restricted 
HF (RHF) calculation. The MO basis allows us to write 
Eq. ^ in a diagonal form equivalent to a more involved 
matrix form when the MO's are developed in an atomic 
basis set [|| of orbitals <p a ,j where a indicates the atomic 
center and j the type: cfo(r) = Y, a ,j C a ,jfa,j(r)- The co- 
efficients Q j, as well as the weights Xj, can be used as 
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variational parameters denning the geminal in Eq. ([6]). 
By truncating the expansion in Eq. ([5]) to a number of 
MO's n equal to the number of electron pairs plus the un- 
paired orbitals, namely n = iVj, one recovers the usual 
RHF theory, because the antisymmctrization operator A 
singles out only one Slater determinant. Moreover, the 
MO weights Xj affect only an overall prefactor of this 
Slater determinant, so that their actual values are irrele- 
vant in this case. However, the pairing function is gener- 
ally not limited to have only N± non vanishing eigenval- 
ues Xj. Therefore, the RVB wave function represents a 
clear extension of the RHF theory, not only for the pres- 
ence of the Jastrow factor, which considerably improves 
the dynamical correlations, but mainly because its de- 
terminantal part goes beyond RHF when n > iV-j-, by 
including also non-dynamical correlations. Quite gener- 
ally, a gain in energy and a more accurate calculation are 
expected whenever n > iVj. 



B. Valence bond energy consistent number of 
molecular orbitals in the AGP 

The main property used in the following derivation re- 
lies on the fact that the atoms are well described by a 
Jastrow correlated RHF (JHF) wave function. Indeed, 
the application of the Jastrow factor J on a simple HF 
Slater determinant provides at least ~ 90% of the correla- 
tion energy in all the atoms (see Refs. @,[T3] and TableU}. 
Here we show that it is possible to extend and remark- 
ably improve the correlated HF approximation for the 
chemical bond, by means of a RVB wave function with 
an appropriate number n of MO's appearing in Eq. (|6j). 
These MO's are chosen to minimize the energy expec- 
tation value in presence of the Jastrow factor, while an 
upper bound on the number n is univocally determined 
by imposing that, when the atoms are at large distance, 
we cannot obtain an energy below the sum of the JHF 
atomic energies. 

The above mentioned criteria are based on the assump- 
tion that the large intra-atomic correlations do not affect 
the chemical properties of the bond, which are instead ex- 
tremely sensitive to the usually much weaker inter-atomic 
correlations. Moreover, electrons close to the atomic cen- 
ters are chemically inert because they are far away from 
the region where the bond is formed. Hence, an im- 
provement in the description of the atoms with many 
determinants [HI would lead in this case only to a rigid 
shift of the total energy. The above assumption is a quite 
generally accepted idea, that has been exploited in differ- 
ent ways by a large variety of approaches. For instance, 
it validates the use of pseudopotentials, the configuration 
interaction (CI) with the frozen core approximation [l^| . 
and is the basis for other quantum chemistry meth- 
ods such as the symmetry-adapted perturbation theory 
(SAPT)01, and the Morokuma analysis (see Ref. H3 
and references therein). 

In the following we shall denote the aforementioned ap- 



propriate number n of MO's with n* . Let us denote with 
M a molecule composed by atoms A\, A2, etc. The opti- 
mal value of n* is most generally obtained by saturating 
a simple upper bound value h: 

r**<n = £iV T (^), (7) 

i 

where i is an index running over the atoms composing 
the molecule M. Since in some cases convergence in the 
energy for the JAGP can be obtained even for n* < h, 
we have used the inequality to define n* and the corre- 
sponding wave function will be denoted by JAGPn*. If 
the sum of the number of spin-up electrons in the atoms 
equals the number of the MO's required by a RHF calcu- 
lation uhf(M) for the molecule, then n* = uhf(M) and 
the JAGPn* wave function reduces to a JHF description 
of the molecule. This is the case, e.g., for Be 2 and B 2 . 
In all the other cases we have uhf(M) < n, and, in this 
work, we have found that there is a substantial energy 
gain in increasing the number of molecular orbitals with 
respect to the RHF value. This happens for instance for 
N 2 , 2 , F 2 , and CN, whereas for LiF, though h > uhf, 
accurate results can be obtained even with n* = tlrf- 

The upper bound in Eq. |(7j) can be slightly improved, 
as it will be shown in the following. This is particu- 
larly important when some degenerate multiplcts of or- 
bitals are not completely occupied, as for the C 2 molecule 
where, by using n molecular orbitals in the AGP ex- 
pansion, one of the two antibonding tt* orbitals remains 
empty, and therefore it is not possible to satisfy the or- 
bital symmetry of the X E+ C 2 wave function. In the gen- 
eral case the highest molecular orbital included in the 
AGP has degeneration D and it may occur that only 
D < D orbitals of the multiplet are included in the AGP 
expansion by the upper bound in Eq. ([7]). For this reason 
it is important to improve the upper bound for n* , in 
particular cases when the chemical compound is spatially 
symmetric, namely for reflections, rotations, translations, 
of the atomic positions. In fact, let us suppose that the 
molecule is composed by several atoms. Some spatial 
symmetry operations can make equivalent ua > 1 iden- 
tical atoms of type A. Assuming that these symmetries 
remain valid up to the atomization limit, we denote by 
m the minimum value of ua among all atomic species. 
Then if m > 1 it is possible to improve the upper bound 
© by: 

n* < n + m- 1. (8) 

For instance for C 2 , according to Eq. ([8]) we have m = 2 
due to the reflection symmetry of the molecule, and 
n* < h + 1. Indeed n* = n + 1 not only allows to fulfill 
the symmetry, but also provides a substantial im- 
provement of the binding with respect to n* = n (see 
Fig. [3] in Appendix IB The one extra molecular or- 
bital added cannot have any effect at large distance in a 
fully symmetric calculation that connects the compound 
at rest to m = 2 equivalent Carbon atoms at large dis- 
tance. Indeed, in this case, the presence of the extra 
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TABLE I: Atomic energies for Li, Be, B, C, N, O, and F: comparison 
between RHF benchmarks, estimated exact values, VMC and LRDMC 
JHF data, and the percentage of recovered correlation energy (%) (eval- 
uated using the estimated exact value and the Hartree-Fock energy). 
For Li and Be all-electron results are shown. For all the other atoms, 
results were obtained with a pseudopotential (Ref. [3). 



Atom 


RHF 


Est. exact 


JHF VMC 


JHF LRDMC 


% 


Li 


-7.432727 a 


-7.47806 a 


-7.47707(6) 


-7.47807(3) 


100% 


Be 


-14.573023 a 


-14.66736 a 


-14.64747(9) 


-14.6575(1) 


89.5% 


B 


-2.54375616 b 


-2.61940948 b 


-2.6031(1) 


-2.6110(1) 


88.9% 


C 


-5.32903005 b 


-5.43249352 b 


-5.4105(1) 


-5.4216(1) 


89.5% 


N 


-9.66837630 b 


-9.79973109 b 


-9.7771(3) 


-9.7898(1) 


92.4% 





-15.70844748 1 


-15.90165954 1 


-15.8754(1) 


-15.89233(8) 


95.2% 


F 


-23.93849161 1 


-24.19290003 1 


-24.1680(3) 


-24.1860(2) 


97.3% 



a From Ref. [TJ 
b From Ref. [l(| 



orbital could improve only the energy of one of the two 
JHF atoms, thus violating their equivalence. Therefore 
the eigenvalue Xj of Eq. ([6]) corresponding to the extra 
molecular orbital must vanish in the atomization limit. 

Generally speaking a value for n larger than the upper 
bound ([8]) certainly leads to a lower value of the total en- 
ergy, but may improve much more the atomic energies, 
rather than the bonding. Actually, we have seen that, in 
all cases so far considered, the accuracy in describing the 
chemical bond improves systematically by increasing the 
number of molecular orbitals, provided it remains smaller 
than the upper bound. Clearly, whenever Eqs. ([7]) or ([5]) 
are satisfied the atomization energy has to be referred 
to the JHF calculation, even if lower energies could be 
achieved with a J AGP wave function for the atoms Q. 
Remarkably, in the limit of large number of molecular or- 
bitals, when the lowest JAGP total energies are obtained 
both for the atom and the molecule, the binding energy 
becomes always worse than the corresponding JAGYn* . 

The JAGPn* wave function can be used also to de- 
scribe bulk systems by applying the upper bound of 
Eq. ([7]) and of Eq. ||5J) to the supcrcell containing a fi- 
nite number of atoms, so that the values of n and m eas- 
ily follow exactly as in the case of a finite open system. 
The upper bound computed in this way may exceed by a 
large amount the number uhf of molecular orbitals nec- 
essary to define a single Slater determinant in the super- 
cell. Thus, convergence in the energy is expected in this 
case for n* <C h. For instance, in the case of graphcnc 
for a typical supercell of 48 atoms, n = 4 X 48 = 192, 
whereas n HF = 3 x 48 = 144 < 192. 

The constrained optimization of the JAGP wave func- 
tion with a given number of MO's is a generalization of 
the standard QMC optimization [l8T| which minimizes the 
total energy, and will be described in Appendix [XJ 



III. RESULTS 

In this section we shall describe the results that we 
have obtained for a set of molecules composed of first row 
atoms, where strong non-dynamical correlation and/or 
weak vdW interactions are present. These molecules are 
used as a test-case for our method. 

Our study has been carried out by means of QMC sim- 
ulations. We started from the constrained optimization 
of the variational wave function described in the previ- 
ous section, which was initialized by taking density func- 
tional theory orbitals in the local density approximation, 
and then we performed variational Monte Carlo (VMC) 
or lattice regularized diffusion Monte Carlo (LRDMC) 
simulations [22| . 

For the determinantal part of the wave function we 
have used a Slater (for Be 2 , and the Li atom in the LiF 
molecule), or mixed Slater/Gaussian (for B 2 , C 2 , N 2 , 2 , 
F 2 , CN, and the F atom in the LiF molecule) basis, large 
enough for an accuracy of lmH in the total energies. 
This quantity sets the tolerance for our complete basis 
set (CBS) limit extrapolation. In particular, for Be 2 the 
basis set is 6s4p2d, for B 2 5s4pld, for C 2 5s5p, for N 2 
5s3p2d, for 2 6s5p2d, for F 2 5s5p2d, for the Li atom in 
the LiF molecule 5s4p, whereas for the F atom, as well 
as for the C and N atoms composing the CN molecule, 
we used the same basis adopted for the corresponding 
dimers. In the mixed Slater/Gaussian cases we have used 
one Slater orbital for each angular momentum, except for 
d orbitals, which have been chosen of a purely Gaussian 
form. Thus, by fully optimizing all the coefficients and 
the exponents of the primitive basis set, we have verified 
that the dimension of the basis is sufficient to achieve 
the desired accuracy. In Appendix IB 21 we show, as an 
example, selected studies of convergence in the basis set. 

A much smaller basis was used for the Jastrow fac- 
tor, because this allows for a more efficient energy opti- 
mization. On the other hand, the essentially exact con- 
tribution of Jastrow-type dynamical correlations, which 
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do not change the phases of the wave function, can 
be very accurately obtained with the well established 
DMC technique|12|. within the recent lattice regularized 
diffusion Monte Carlo (LRDMC) implementation^. 
LRDMC is equivalent to standard DMC for all-electron 
calculations, and represents an improvement of the older 
technique because it allows to obtain a rigorous upper 
bound of the total energy even when pseudopotentials 
are used in the calculation. The DMC /LRDMC ap- 
proach can be seen as a stochastic optimization of the 
Jastrow factor which keeps fixed the phases of the RVB 
wave function. In some test cases (see Appendix IB 2[) , 
we have also verified that a larger basis in the Jastrow 
does not provide significant changes in the physical and 
chemical quantities here considered, because total energy 
differences are much less sensitive to the extension of the 
Jastrow basis set. 

We have used a helium-core pseudopotential[l6| for all 
but Be and Li atoms. In some test cases without pseu- 
dopotentials (e.g. Be 2 ) we have explicitly verified that 
the DMC and the LRDMC energies are consistent, but 
we have adopted the latter method for the sake of gener- 
ality. In the C 2 case we have also checked that the effects 
of the pseudopotential on the total energy differences are 
negligible [23j. 

In Table UH we compare with estimated exact results 
bond lengths and well depths obtained by means of 
VMC and LRDMC simulations performed with JHF or 
JAGPn* wave functions for the various molecules consid- 
ered in this paper. We optimized each wave function for a 
bunch of different interatomic distances. The energy and 
interatomic distance at the minimum were found by in- 
terpolating the energy close to its minimum value with a 
cubic polynomial. We also report binding energies found 
in Ref. [18| with DMC calculations for a fully optimized 
all-electron Jastrow-correlated single determinant wave 
function (in our table denoted with JxSD DMC). Finally, 
we compare the JAGPn* zero point energy (ZPE) with 
available experimental data. This quantity was com- 
puted by standard first order perturbation theory in the 
anharmonic cubic term. For this property, the agreement 
between both VMC and LRDMC results and experimen- 
tal findings is satisfactory in most of the cases. The ac- 
curacy in the ZPE can be probably improved by doing a 
more careful fit around the minimum. 

Below we comment all the different cases. 



A. Beryllium and boron dimers 

Though the Beryllium dimer does not belong to the 
so-called 'Gl-set' of molecules (see e.g. Ref. 14]), this 
dimer is a very important test case for the variational 
method proposed in this paper. Indeed, several compu- 
tational methods (see e.g. Ref. [3(|), including previous 
QMC simulations Q, have failed in the attempt of re- 
producing the binding of this molecule. Moreover, until 
the '80s Be 2 represented a technical challenge for the ex- 



perimentalists, and even later the value of its binding 
energy was not well established. A review on the experi- 
mental and theoretical investigations of Be 2 has recently 
appeared [24|], containing also new reference experimental 
data for its binding energy. 

In Fig. [T] we provide the energy dispersion curve for 
the Be 2 molecule. The main plot shows a compari- 
son between standard RHF calculations [3(J, VMC data 
obtained with the JAGPn* wave function, VMC and 
LRDMC results for a J AGP with n > n* . We also show 
an expanded Morse oscillator (EMO) fit of the recent 
experimental data of Ref. [24[ . 

As mentioned before, in this case it turns out that 
n* < n#i?(Be 2 ) + 1. In particular, by using the upper 
bound of Eq. Q our JAGPn* reduces to a simple JHF 
wave function with n* = 4 [the upper bound n* = 5 
of Eq. ([8]) does not provide significant improvements 
in a fully symmetric calculation]. Within the n = n* 
constraint, bond features such as binding energy and 
bond length are reproduced fairly accurately, whereas 
a trial wave function with n > n#i?(Be 2 ) + 1 fails to 
bind the molecule at the expected distance, even though 
the total VMC (LRDMC) energy E = -29.32295(8)H 
(E = — 29.33385(7)H) is much below the constrained 
minimization by about 24mH (14mH) at R = 5 a.u.. 
This total energy is very accurate from an absolute point 
of view and compares well with state of the art QMC 
calculations. [18[ However the variational wave function 
with the lowest variational energy, i.e. the JAGP with 
n = 10, behaves similarly to an uncorrelated RHF, and 
both provide a very poor description of this chemical 
bond.[3l| More in detail, the VMC JAGPn = 10 en- 
ergy dispersion curve presents one minimum at an inter- 
atomic distance R > 8 a.u., while LRDMC JAGPn = 10 
displays an additional swallower minimum close to the 
expected bond length. On the other hand, the quite ac- 
curate dispersion curve obtained by the full optimiza- 
tion of the JHF wave function shows, for the first time 
to our knowledge, that the key missing ingredients in 
the HF for Be 2 are just the dynamical correlations car- 
ried out by our Jastrow factor. Though very simple, 
our Jastrow factor includes many-body correlations (up 
to two- ion two-electron interactions), that allow to take 
into account effective attractions between atoms given 
by vdW forces, 13] and other polarization-polarization 
contributions. [32| Indeed, the dynamical interactions are 
extremely important to bind the molecule and it is crucial 
that the Jastrow factor includes this effect. For instance, 
the different parametrization of the Jastrow factor used 
in Ref. [18| does not allow to bind Be 2 at a variational 
level, at variance with this work. On the other hand, 
the DMC binding energies of Ref. [l8[ are much closer 
to our VMC and DMC results, further suggesting the 
importance of the dynamical correlations in the bond. 

In the inset of Fig. [T] we compare the VMC and 
LRDMC JAGPn* energy dispersion curves shifted by 
their asymptotic limits. Despite some slight differences, 
the agreement between the two QMC techniques within 
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TABLE II: Bond lengths (in a.u.), well depths (in eV) and ZPE (in mH) for a set of first row 
diatomic molecules. We report VMC and LRDMC values for both JHF and JAGPn* trial wave 
functions, and experimental results or estimated exact values. The well-depth exact estimates are 
given by the experimental binding energies subtracted by the spin-orbit energies when accessible 
(i.e. for all atoms but B 2 and Be 2 ) and the ZPE. We also report the JxSD DMC well depths 
of Ref. pl| when available. For Be 2 (all electron calculations) and B 2 (calculations with the 
pseudopotential in Ref. [HI]), n* = JVf = 4, hence the JHF and JAGPn* results coincide. 







B 


9nd length (a.u.) 












Be 2 (all el.) 


B 2 


c 2 


N 2 


o 2 


F 2 


LiF 


CN 


JHF VMC 


4.85(5) 


3.041(6) 


2.367(2) 


2.048(1) 


2.27(1) 


2.66(1) 


2.95(4) 


2.185(6) 


JAGPn* VMC 


4.85(5) 


3.041(6) 


2.334(6) 


2.075(2) 


2.268(7) 


2.661(5) 


2.92(2) 


2.200(6) 


JHF LRDMC 


4.65(7) 


3.021(9) 


2.369(3) 


2.051(1) 


2.270(4) 


2.665(9) 


2.949(8) 


2.201(3) 


JAGPn* LRDMC 


4.65(7) 


3.021(9) 


2.337(6) 


2.075(1) 


2.277(4) 


2.663(3) 


2.950(7) 


2.202(2) 


Exact estim. 


4.63 a 


3.005 b 


2.3481 c 


2.075 b 


2.283 b 


2.668 b 


2.955 b 


2.214 b 


Well depth (eV) 




Be 2 (all. el) 


B 2 


c 2 


N 2 


o 2 


F 2 


LiF 


CN 


JHF VMC 


0.120(5) 


2.754(3) 


5.538(9) 


9.662(3) 


4.976(8) 


1.124(4) 


5.93(2) 


7.52(1) 


JAGPn* VMC 


0.120(5) 


2.754(3) 


6.327(9) 


9.874(2) 


5.060(7) 


1.671(2) 


5.96(2) 


7.68(1) 


JxSD DMC 


0.125(1) 


2.798(3) 


5.656(3) 


9.583(3) 


4.992(7) 


1.349(6) 






JHF LRDMC 


0.143(6) 


2.797(2) 


5.763(9) 


9.665(2) 


5.070(5) 


1.452(3) 


6.049(6) 


7.661(5) 


JAGPn* LRDMC 


0.143(6) 


2.797(2) 


6.297(8) 


9.882(1) 


5.126(5) 


1.688(2) 


6.056(6) 


7.744(5) 


Exact estim. 


0.1153(3) a 


2.91(6) d 


6.43(2)° 


9.902(3) c 


5.233(3)' 


1.693(5) c 


6.03(9) f 


7.86(9) f 


ZPE (mH) 




Be 2 (all el.) 


B 2 


c 2 


N 2 


o 2 


F 2 


LiF 


CN 


JHF VMC 


0.56(5) 


2.49(5) 


4.3(1) 


6.38(6) 


3.8(1) 


2.20(3) 


2.3(2) 


4.9(1) 


JAGPn* VMC 


0.56(5) 


2.49(5) 


4.2(1) 


5.48(3) 


3.85(9) 


2.20(3) 


2.1(2) 


4.87(8) 


JHF LRDMC 


0.61(9) 


2.51(7) 


4.38(3) 


5.83(6) 


3.77(5) 


2.16(3) 


2.18(8) 


4.81(3) 


JAGPn* LRDMC 


0.61(9) 


2.51(7) 


4.3(1) 


5.51(2) 


3.70(9) 


2.22(2) 


2.10(6) 


4.82(4) 


Exp. 


0.56 a 


2.4 b 


4.2 C 


5.4 e 


3.6° 


2.1 e 


2.07 b 


4.71 b 



a From Ref. O 

b From Ref. [H 

c From Ref. [2(| 

d From Ref. Hll 

c From Ref. 1281 

From Ref. [H 



the n = n* constraint and the most recent experimental 
fmdmgs[24| can be considered fairly good in this case, 
due to the very weak binding of the molecule. 

Also the JAGPn* description for B 2 reduces to a JHF 
wave function. Both bond length and binding energy 
agree within two standard deviations with the estimated 
exact data. 



B. Fluorine dimer 

A remarkable example of the accuracy of our technique 
is provided by the energy dispersion curve of the fluo- 
rine dimer reported in Fig. [21 where we show the re- 
sults obtained with various QMC methods (and differ- 
ent wave functions), and other ab-initio results. More 
in detail, we compare our JHF and JAGPn* VMC data 
(see also Table IIII1 where, for comparison, we also re- 
port our LRDMC results) with two energy dispersion 
curves obtained with auxiliary-field QMC (AFQMC) 
simulations for an unrestricted HF reference wave func- 
tion spin-projected to eliminate spin contamination [33| , 
and an ab-initio study based on full configuration in- 



teraction (FCI) calculations combined with the correla- 
tion ener gy extrapolation by intrinsic scaling (CEEIS) 
technique [3J] plus core-electron correlations and scalar 
relativistic corrections [35j. 

One can observe the dramatic improvements of the 
JAGPn* wave function with respect to JHF simulations 
(see also Table HH - Tni[) . According to Eq. ([7]), we have 
used n* — n^fF(F 2 ) + 1, because the upper bound of 
Eq. ©, n* = nffF(F 2 ) + 2, does not lead to signifi- 
cant differences within our energy accuracy. We remark 
here instead the importance of adding just one molecu- 
lar orbital to the Hartree-Fock theory, because this al- 
lows to consider all bonding and antibonding MO's in 
the AGP, thus leading to a fully size consistent result 
which benchmarks the energy dispersion curve from the 
bond length to the atomization limit. The agreement of 
the JAGPn* with the ab-initio CEEIS-FCI calculations 
is remarkably good already at a VMC level. In fact, 
the VMC binding is 1.671(2) eV against 1.6867 eV of 
FCI calculations (without spin-orbit corrections). The 
LRDMC binding is 1.688(2) eV. Instead the AFQMC 
curves seems to be shifted of approximately 2-3mH with 
respect to our JAGPn* data in the bond and intermedi- 
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TABLE III: F 2 binding energies (in mH) as a function of the internuclear 
distance R in a.u. (see also Fig. [5J): VMC and LRDMC results (energy 
of the molecule at distance R minus two times the JHF atomic energy) 
for a JHF and the JAGPn* wave functions. 



R (a u 1 


VMC JHF 


VMC JAGPn* 


LRDMC JHF 


LRDMC JAGPn* 


2.36 


-23.9(4) 


-39.9(4) 


-32.4(4) 


-39.6(4) 


2.46 


-28.1(6) 


-52 3f5^ 


-42.1(5) 


-53.2(4) 


2.56 


-35 n(nS 


-59 


-48.2(5) 


-59.7(4) 


2.668 


-39 6(5^ 


-61.0(4) 


-49 l(i') 


-60.9(4) 


2.76 


-35.6(5) 


-59.7(5) 


-49 1(4) 


-60.6(4) 


2.86 


-30.1(4) 


-56.2(4) 


-44.1(4) 


-57.1(4) 


2.96 


-22.0(5) 


-51.9(4) 


-38.1(5) 


-52.1(4) 


3.3 


3.9(5) 


-34.2(5) 




-34.6(4) 


3.8 


36.5(6) 


-14.4(5) 




-13.4(4) 


4.5 


41.2(6) 


-4.0(4) 




-2.5(4) 


5.5 


28.1(5) 


-1.0(4) 






6.5 


12.8(5) 


-0.1(5) 






7.5 


8.5(5) 


0.5(5) 






8.5 


5.4(5) 


0.2(5) 







ate length regions. This is due to an underestimation of 
the energy at large distance caused by the use of a simple 
unrestricted HF wave function (see the discussion in Sec- 
tion IV of Ref. [HI and Figs. 4 and 6 therein). Indeed the 
AFQMC well depth is 1.70(2) eV and 1.77(1) eV for the 
cc-pVTZ and the cc-pVQZ wave functions respectively, 
when the reference at large distance is the molecular en- 
ergy, whereas it is 1.60(1) eV and 1.70(1) eV respectively, 
when the large distance reference is twice the energy of 
the separated atoms. 

The results we have presented so far, reported in Ta- 
ble HH represent an astonishing example of the impor- 
tance of constraining the variational wave function to an 
appropriate form during the optimization of the energy. 
Indeed, a brute force optimization of a correlated wave 
function, which is a rather demanding computation espe- 
cially within QMC, would lead to an upper bound of the 
total energy which is almost meaningless, particularly in 
the Be 2 case. The rational behind this effect is that an 
unconstrained optimization may not satisfy the require- 
ment for the wave function to be a fair representation of 
the ground state of a physical Hamiltonian. While in lat- 
tice models it is possible to constrain the determinantal 
part of the RVB wave function to be the ground state of 
a short-range Bardeen-Cooper-Schrieffer (BCS) Hamilto- 
nian -a quite sensible and accepted choice in strongly cor- 
related models- this is much harder in continuous-type 
calculations. The constraint that we propose, very simple 
to implement in practice, just mimics the effect of com- 
puting the ground state of an HF Hamiltonian with an 
additional sufficiently weak BCS coupling between elec- 
trons. In fact, in this limit one obtains the complete or 
partial occupation -via the Xj in Eq. of a number of 
molecular orbitals n* not necessarily equal to the RHF 
prediction. In this context, the BCS coupling represents 
the effective interaction between electrons, which pairs 
them into the chemical bond. For instance, it is well 
known that the ground state of the H 2 dimer at large 



distance is very well described by the singlet entangled 
state obtained with the AGP[36|. only when the bond- 
ing and antibonding orbitals are taken into account. This 
state can be considered the ground state of a BCS Hamil- 
tonian that in the atomization limit simply splits into a 
sum of two atomic HF Hamiltonians, with vanishingly 
small pairing. This coupling is however important to lift 
the degeneracy between the singlet and the triplet states. 
The same physics happens in the F 2 molecule studied in 
this work. 

C. Carbon, nitrogen, and oxygen dimers 

The C 2 , N 2 and 2 molecules represent challenging 
cases for our correlation consistent AGP approach. In- 
deed, when < S < S max and S max > 1, the lack of 
size-consistency in the AGP poses a fundamental limita- 
tion in order to reach the JHF limit in the dissociation. 
Note that S = and S max = 1 is a non trivial case when 
the JAGP is size consistent and the JHF is not (e.g. the 
simplest H 2 molecule, or the F 2 described in the previous 
section). Strictly speaking, the restriction of the number 
of MOs to n* does not guarantee a size consistent JHF 
result, even for the JAGP wave function. In the general 
case, the total JAGP energy in the atomization limit is 
an upper bound of J2a Ejhf{A) evaluated in the CBS 
limit. A generalization of the JAGP, based on the pfaffian 
algebra, which includes also triplet pairing for electrons 
with the same spin, allows to have JHF size consistent re- 
sults also in the cases with S = S max — 1 (S max > 1), e.g. 
in 2 . Despite we have not implemented this generaliza- 
tion, triplet pairing seems to provide a rather negligible 
effect at bond length, as very good results can already be 
obtained with the present JAGP ansatz. 

The constrained minimization of the JAGPn* wave 
function leads to very significant improvements with re- 
spect to JHF results in both the binding energy and the 
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4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 
R (a.u.) 



FIG. 1: Be 2 binding energy (in milli-Hartree) versus the nu- 
clei distance R (in atomic units). Comparison of RHF[30| 
outcomes, JAGPn* results (in this case n* =4), and VMC 
and LRDMC results obtained with a JAGP wave function 
with n — 10 (squares, dots, downward and upward triangles 
respectively, while lines are a guide to the eye). In the figure, 
the experimental binding energy subtracted by the zero point 
energy [241] (solid line), and an EMO fit of the experimental 
data[2J] (slash-dotted line) are also plotted. The reference 
atomization limit for the JAGPn* results is given by atomic 
calculations with a JHF wave function. For n — 10, the atom- 
ization reference is given by an atomic JAGP wave function 
with the same primitive basis set as the JAGP wave func- 
tion for the dimer. In the inset: comparison of fits to VMC 
and LRDMC data (solid and dashed lines respectively) for 
the JAGPn* wave function and the EMO fit of experimental 
data from Ref. [24J (slash-dotted line). Labels for the inset 
axes are the same as in the main frame. All curves are shifted 
with respect to their own atomization limit. 







Exact estim. 








VMC JHF 








CEEIS FCI 


A 






VMC JAGPn* 






V" I 


AFQMC cc-pVTZ 


— v — 




AFQMC cc-pVQZ 













2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 
R (a.u.) 



FIG. 2: Energy dispersion curve for the F 2 molecule obtained 
with several computational techniques. The zero reference en- 
ergy (dashed line) is twice the JHF atomic total energy for the 
JAGPn* and the JHF molecules, whereas it is the large dis- 
tance energy of the dimer for the remaining data. The slash- 
dotted line indicates the experimental binding energy[2g] sub- 
tracted by the zero point and spin orbit energies. Lines are 
a guide to the eye. The CEEIS-FCI data are taken from 
Ref. |H, whereas AFQMC data are taken from Ref. EH- 



TABLE IV: Well depths for C,, N 2 2 : 
comparison between pfaffian results [loL KL1I] and 
JAGPn* results. The VMC and LRDMC find- 
ings are compared with the exact estimates pre- 
viously reported in Table [IT] 



Method WF 


c 2 


N 2 


o 2 


VMC STU 
VMC JAGPn* 
DMC STU 
LRDMC JAGPn* 


5.94(2) 
6.327(9) 
6.26(2) 
6.297(8) 


9.42(3) 
9.874(2) 
9.84(2) 
9.882(1) 


4.94(3) 
5.060(7) 
4.93(2) 
5.125(5) 


Exact, estim. 


6.44(2) 


9.908(3) 


5.241(3) 



bond length for C 2 and N 2 , if compared with the exact es- 
timates coming from the experimental values. In 2 the 
situation is quite different, since n* = nup(0 2 ) + 1 and 
the molecular orbital missing in the HF scheme is quite 
high in energy, compared to all the paired molecular or- 
bitals included in Eq. ([(TJ according to our constraint. 
Therefore, it is not surprising that, in this case, the im- 
provement upon JHF findings is smaller. The comparison 
with the exact estimates is nevertheless quite satisfactory. 

By comparing the results in Table [TT1 for N 2 , 2 and 
F 2 with Ref. [l^j's single-determinant DMC well depths, 
we note improvements already at the JHF level, and 
even bigger improvements are obtained upon Grossman's 
benchmarks [1J] (to which one should add the ZPE). We 
should mention that in Ref. [3, the basis used for the 
determinantal part is not optimized. Other differences 
could be due to either the different pseudopotentials used 
or the calculations not fully converged in the CBS limit. 

As a further evidence of the accuracy of our method, 
in Table [IV] we compare the well depths found for C 2 , 
N 2 and 2 with results obtained by using a singlet- 
triplet-unpaired (STU) pfaffian wave function[l(| HH 
and taking as a reference at large distance the JHF 
atomic limit, i.e. the binding energy is computed as 
E molecule {ST\J) -2E atom (mF). As shown in Table El 
we find a rather good agreement for the well depth even 
in the challenging C 2 molecule, when S — S max — 2 and 
an exact size consistent result is not possible even within 
the more general pairing function containing triplet cor- 
relations. In principle, the pfaffian wave function, being a 
generalization of the JAGP, should have a larger binding 
energy by taking the corresponding JHF atomic energy 
as a reference, provided the CBS limit is reached and the 
pseudopotentials are accurate enough. Instead, in all the 
cases shown in Table \N\ the STU and the JAGPn* bind- 
ing energy are very close, at least at the DMC level. The 
fact that our binding energy is always larger and more 
accurate comes probably from the use of a more com- 
plete basis, fully optimized in both the coefficients and 
the exponents, whereas in Ref.[l(| EH the atomic basis 
is not optimized. 

These comparisons show that our JAGPn* ansatz pro- 
vides quite generally a very accurate description of the 
chemical bond. This emphasizes the role of singlet elec- 
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tron pairs in the chemical bond, and is consistent with the 
RVB theory. Thus, the present variational wave function 
can be considered the cheap but nevertheless accurate re- 
alization of the RVB idea, since only a single determinant 
and standard variational Monte Carlo are needed. 



D. Heteronuclear dimers 

We further carried out calculations for a couple of het- 
eronuclear diatomic molecules belonging to the Gl set 
(namely LiF and CN), selected on the basis of the quite 
big discrepancy between the binding energy found in 
Ref. [lH and the reported experimental values. Bond 
lengths, dissociation energies and ZPEs for these two 
molecules are reported in Table [Til We compare our well 
depths with the exact estimates obtained by correcting 
the experimental dissociation energies with the experi- 
mental ZPE and spin orbit energies reported in Ref. [29l | . 

As mentioned in Section Hi Bl LiF is one of the cases in 
which n* — n does not yield significant improvements 
with respect to the JHF wave function, even though 
h > uhf- Instead, for CN the JAGPn* wave function 
improves the description of the bond with respect to the 
JHF one, giving a bond length in fairly good agreement 
with the experimental value, although the binding en- 
ergy is underestimated by « 0.1 eV with respect to the 
exact estimate reported in Ref. [2{J. As for the homonu- 
clear dimers shown in the previous sections, also in the 
heteronuclear cases here considered, our method provides 
binding e nerg ies in closer agreement with the experimen- 
tal values [291] than those of Ref. [l4| . In particular, for 
LiF the agreement is very good already at the JHF level, 
as anticipated above. 



IV. CONCLUSIONS 

In conventional quantum Monte Carlo variational tech- 
niques, based on the use of the Jastrow factor, it is not 
possible to consider a finite basis set and to exploit the 
huge cancellation between atomic energies and molecu- 
lar energies within the same basis set. Indeed, after the 
introduction of the Jastrow factor, the wave function is 
unavoidably defined on an infinite dimensional Hilbcrt 
space. As a consequence, it is more difficult to achieve 
the chemical accuracy on the energy differences and ob- 
tain a good description of the chemical bond, as we have 
shown for instance in the Be 2 case. Here, a very accurate 
variational energy obtained by applying the DMC tech- 
nique to our lowest energy JAGP wave function, com- 
pletely misses the correct features of the bond. In this 
case, with an unconstrained variational approach, quali- 
tatively correct results can be obtained probably only by 
reaching the chemical accuracy on the total energy, that 
is clearly a very difficult task for any approximate varia- 
tional technique. In fact, this target was so far achieved 



within QMC only by using several determinants in small 
molecules [18{. 

In this paper we propose a simple constraint which al- 
lows to exploit the above mentioned cancellation between 
atomic and molecular energies even in QMC calculations 
based on a single determinant wave function. In fact, 
instead of imposing a constraint on the dimension of the 
atomic Hilbert space we change a bit this point of view 
by constraining the number of molecular orbitals to an 
appropriate value that allows to take JHF results for the 
isolated fragments as a reference for the dissociation en- 
ergy of the molecule. With this constraint we have shown 
that it is possible to obtain much more accurate results 
in both variational and LRDMC calculations. 

Although we have not carried out a systematic study of 
all the Gl set considered in Ref. [3 in several cases where 
the discrepancy was sizeable we obtain an almost exact 
description of the bond (e.g. in F 2 ). Surprisingly the 
LRDMC calculation provides only small improvements 
upon the simple and much cheaper VMC calculation, 
which turns out to be remarkably accurate in our ap- 
proach. Also in cases where we do not improve upon 
the JHF results (e.g. LiF), we nonetheless obtain accu- 
rate binding energies, within a precision of about 0.1 eV. 
The latter achievement could be due to the accurate ba- 
sis set we have considered in our work, together with the 
state-of-the-art optimization technique. [371] which is able 
to handle the large number of parameters in an extended 
basis set. 

In conclusion, we have introduced a new and general 
approach to perform electronic structure calculations of 
quantum chemistry compounds based on a variational 
RVB wave function. In this formulation, we have shown 
that a substantial improvement in the description of the 
chemical bond is possible by extending the standard cor- 
related single determinant theory with the JAGP wave 
function. In the original formulation of the RVB theory, 
the gain in energy obtained by the resonance of several 
valence bond configurations was just named the 'reso- 
nance valence bond energy'. Within this new formula- 
tion we propose that this energy gain can be achieved by 
increasing the number of molecular orbitals of the JAGP 
from its HF value, and without exceeding a value n* of 
molecular orbitals. This value can be determined by re- 
quiring a correlation consistent property from the bond 
length to the atomization limit, realized via a constrained 
energy minimization. 
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APPENDIX A: CONSTRAINED OPTIMIZATION 
OF THE AGP WAVE FUNCTION 

1. Molecular orbital expansion of the AGP 

In this appendix we expand the pairing function $ in 
atomic orbitals £j (r) located at atomic positions Rj : 



(Al) 



where A is the pairing matrix and j, f label the consid- 
ered atomic orbitals on the corresponding atomic posi- 
tions Rj , Rj/ . Obviously, in order to define a singlet state 
the pairing matrix should be symmetric Aj.j' = Xj'j- 
Hereafter, both for simplicity and for the sake of gen- 
erality we do not assume this symmetry, because it can 
be easily satisfied during the optimization scheme, when 
necessary. Therefore, in the general case we are left with 
Nl = L x L — 1 independent variational constants, where 
L is the linear size of the matrix A, namely the dimen- 
sion of the atomic basis. There is only one linear depen- 
dence between the L x L entries of the matrix A because 
the multiplication of by an overall constant does not 
change the AGP apart for its normalization. This con- 
straint is usually satisfied by keeping fixed an arbitrary 
matrix element to the unit value. 

Usually, the number Nl is very large and in the follow- 
ing we determine a systematic way to work with much 
less variational parameters, being nevertheless efficient in 
determining the lowest energy molecular orbitals of the 
chosen variational ansatz. 

For simplicity we do not consider unpaired orbitals, 
because for them no constraint is applied, therefore we 
set iVf — N±. Moreover, in the following we can assume 
that the original orbitals £j have been orthogonalized by a 
suitable transformation depending on the overlap matrix 



m 3 



(A2) 



namely we implicitly assume the following change of the 
definition of the orbitals and the corresponding matrix A 
in Eq. (|AT|) : 



k,l 



(A3) 



This greatly simplifies the forthcoming analysis without 
loss of generality. 

Then, for the resulting square matrix we can use the 
well known singular value decomposition: 



\j — 



fc=i 



(A4) 



where ctk > and -0 fc (i/> fc ) are a set of r < L molecular 
orbitals for the spin-up (spin-down) electrons that are 



orthonormal, i.e. X)z V'/V'/ = S% j. Formally the spin- 
up molecular orbitals and the spin-down ones are the 
eigenvectors of the 2L x 2L symmetric matrix 



H = 



A 
At 



(A5) 



which has pair of eigenvectors with eigenvalues ±y / Ofc 
given by: 



(A6) 



A simple way to reduce the number of parameters is 
to require that the matrix has rank r < L so that all 
the eigenvalues otk for fc > r are assumed to be zero or 
negligible. For instance if r — N/2 we obtain the stan- 
dard Slater determinant with N-\ = N± = N/2 molecular 
orbitals for each spin component. 

This projection scheme can be made general, and this 
leads to a remarkable extension of the Slater determinant, 
within the AGP wavefunction expanded in molecular or- 
bitals, as discussed in the forthcoming subsection. 



2. Projection on a rank-r geminal 

If the rank r of a geminal matrix Ajj is equal to half 
the number of electrons N/2, then the AGP represents a 
Slater determinant. Even if N/2 is usually much smaller 
than the dimension of the atomic basis L, Fermi statistics 
at zero temperature favors the occupation of the lowest 
possible energy levels, so that r ~ N/2 turns out to be 
a reasonably accurate guess for the AGP wave function. 
In principle this wave function may have much larger 
rank up to r = L, but one may expect that most of the 
singular values will have negligible weight. Therefore, 
from a general point of view, and not only for reducing 
the number of variational parameters, it is important to 
optimize in an efficient way a full L x L matrix of rank- 
N/2 < r < L given by Eq. l[A4]) . 

To this purpose we propose the following scheme of 
constrained optimization, where r is chosen and fixed to 
a reasonable value n* « N/2 during the optimization. 

Given A a rank-r matrix, in order to simplify the nota- 
tions, we write the corresponding singular value decom- 
position (|A4[) in a matrix form: 



A = V°a°V4, 



(A7) 



where ip and ip are L x r matrices, the subscript T 
indicates the transpose of a matrix, and the non zero 
singular values a° , k = 1 , ■ • • r are denoted by a diagonal 



matrix a . 

Then we change this matrix A by adding to it a general 
first order contribution: 

A' = A + A 1 (e), 

where henceforth the superscript indicates the order of 
the expansion in e. This new matrix will be constrained 
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to have rank r. Therefore all the terms in Eq. (|A4[) can 
be expanded within first order in perturbation: 

A 1 = £(fa°^ + ^aVT) 

+ e^a 1 ^. (A8) 

In order to satisfy the constraint on the rank in the ma- 
trix A , it is much simpler to work with an unconstrained 
matrix A, and left and right projection matrices: 

P R = (A9) 
P L = (A10) 

The two matrices above are projection matrices (P = Pr 
and P 2 = P) as they project vectors in the r dimensional 
subspaces corresponding to the non zero values of the 
singular value decomposition (|A4[) . 

Indeed it is very simple to show that if the matrix 
A + A 1 satisfies the constraint of a singular value decom- 
position with rank r, A 1 has to satisfy the simple relation: 

(I- P L )X 1 (I-P R ) = 0, (All) 

because in the expression (|A8|) (I — P L )ip° = and 
^°(7 - P R ) = 0. 

Thus an unconstrained variation of the matrix A can 
be projected onto the constrained one by using the above 
projection matrices: 

A 1 = A - (I - P L )X(I - P R ) (A12) 

in the sense that, after the above projection, the matrix 
A° + A 1 is suitable and can be considered to satisfy the 
constraint of a rank-r matrix at first order in the pertur- 
bation (the matrix A 1 being sufficiently small). 

Indeed, by simple inspection, the RHS of Eq. (|A12[) 
immediately satisfies the condition (|A11|) . that is so far 
considered a necessary condition. It is also possible to 
show with a lengthy but straightforward calculation using 
first order perturbation theory of the symmetric matrix 
H given in (|A5|) . that relation (|A1 1|) is also a sufficient 
condition for a perturbation that does not change the 
rank of a singular value decomposition. 

3. Application to QMC 

In the actual application of the recent QMC scheme 
for minimization of the energy, it is important to evalu- 
ate derivatives of a function with respect to the uncon- 
strained parameters A. This function E(x) can be cither 
the logarithm of the wave function or the local energy 
on a particular electronic configuration x sampled by the 
MC technique. @, i,[l| . 

Given the matrix Dij = dE/dXi^ of the uncon- 
strained derivatives with respect to Xij, by using 
Eq. (|A12|) and the chain rule for derivatives, then the 
corresponding matrix of constrained derivatives Dij — 
dE/dXij can be computed by simple matrix manipula- 
tion in the following way: 



D = D - (I - Pt)D(I - Pfi). (A13) 

In order to work with the original matrices we have to 
replace in Eqs. (|A12|A13|) the ones obtained by applying 
the inverse of the transformation (|A3j) : 

pR ^ S l/2pR s -l/2^ (AU) 
pL ^ s -l/2pL 3 l/2_ ( M5 ) 

Notice also that after this transformation P R and P^ are 
no longer equal to P R and P L in Eq. (|A13[) . 

The scheme therefore can be summarized in the fol- 
lowing steps: 

1. Compute the unconstrained derivatives D that, 
with some algebra, can be casted into a product of 
much smaller rectangular matrices U, V such that 
D = UtV of dimension L x N/2. Notice that also 
the projection matrices can be written in this con- 
venient form, as in Eq. (IA9I) . 

2. Apply the projection (|A13|I by using the current 
molecular orbitals. By exploiting the fact that 
all the matrices involved are written in terms of 
much smaller rectangular matrices, a very conve- 
nient computation can be achieved scaling like N 2 L 
instead of L 3 as in the straightforward implemen- 
tation of the projection. 

3. Apply the recent optimization schemes [1, fl8| . 
and change the unconstrained parameters A. Then 
apply the projection, by diagonalizing the matrix 
A and taking only the right and left eigenvectors 
corresponding to the largest singular values. New 
molecular orbitals are then defined after this diag- 
onalization. 

4. Repeat the above-described steps until convergence 
in the energy is achieved. 

APPENDIX B: AGP AND BASIS SET 
EXPANSION EFFECTS 

1. Effect of the improved upper bound for n*: the 

C 2 case 

As explained in Section lll Bi C 2 is one of the exceptions 
to the rule of Eq. ([7]). In this case, we have n = 6, but the 
more accurate upper bound in Eq. ([5]) allows to work with 
n* = 7. Indeed, by following strictly Eq. (|7|) one would 
include a single antibonding orbital 7r* in the AGP, while 
that orbital is double degenerate, due to the rotational 
symmetry of the molecule. Therefore, in order to fulfill 
the symmetry of the dimer, it is particularly important to 
fill the degenerate levels in the AGP by setting n* = 7. 
In Fig. [3] we show the VMC and LRDMC energies at 
various internuclear distances found with a JAGP wave 
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FIG. 3: Energy dispersion curve for C 2 dimer with n = 6 
(circles), n = 7 (squares), n = 8 (triangles) molecular orbitals. 
VMC and LRDMC results are represented by filled and empty 
symbols respectively. 



TABLE V: Total energies from all electron calculations for 
C 2 with a JAGP expanded on a different number of MOs 
n are shown, n — 9 corresponds to the n* upper bound of 
Eq. JS). The estimated exact total energy Eq is reported for 
comparison. 



n 


VMC 


LRDMC 


9 


-75.8439(7) 


-75.8934(4) 


10 


-75.8453(7) 


-75.8930(4) 


15 


-75.8473(8) 


-75.8928(4) 


Eq 


-75.9265 a 



From Ref. 28 



function expanded in h and n* molecular orbitals. In 
this case, the improved upper bound for n* yields a gain 
of w 2.7 mH in the VMC energies and of w 1.5 mH in 
the LRDMC ones. Incidentally, the n = n* = 7 energies 
agree within the error bars with the data resulting from 
a JAGP with n = 8 molecular orbitals. 

An analogous check was done with all electron sim- 
ulations at a fixed interatomic distance R = 2.35 a.u.. 
Results are reported in Table IVl We note a saturation of 
LRDMC total energies for n> n*. 



2. Convergence in the basis set of the AGP and 
the Jastrow parts 

Below we report the convergence in the basis set for 
selected molecules. 



a. N s convergence in the basis set 

For N 2 we checked the convergence in the basis set by 
means of VMC and LRDMC simulations at the experi- 



TABLE VI: JHF and JAGPn* VMC and LRDMC total ener- 
gies (in Hartree) for the N 2 pseudo-molecule. Energies were 
computed at the experimental bond length using different ba- 
sis sets for both the determinantal part and the three and 
four-body Jastrow term (34BJ). 



JHF 


Det. 34BJ 


VMC LRDMC 


5s5p 2s2p 
5s5p 3s3p 
5s5p 4s4p 
4s4p 3s3p 
6s6p 3s3p 
5s3p2d 2s2p 
5s5p2d 3s3p 


-19.9031(5) -19.9338(2) 
-19.9071(3) -19.9346(2) 
-19.9076(4) -19.9349(2) 
-19.9071(3) -19.9344(2) 
-19.9080(3) -19.9348(2) 
-19.9074(3) -19.9339(1) 
-19.9086(4) -19.9349(2) 


JAGPn* 


Det. 34BJ 


VMC LRDMC 


5s5p 2s2p 
5s5p 3s3p 
5s5p 4s4p 
4s4p 3s3p 
6s6p 3s3p 
5s3p2d 2s2p 
5s5p2d 3s3p 


-19.9159(5) -19.9422(2) 
-19.9204(3) -19.9433(2) 
-19.9200(3) -19.9423(2) 
-19.9185(3) -19.9418(2) 
-19.9205(3) -19.9430(2) 
-19.9208(3) -19.9430(2) 
-19.9213(3) -19.9430(2) 



TABLE VII: JHF and JAGPn= 8 total energies (in Hartree) 
for the C 2 pseudo-molecule. Energies were computed at the 
experimental bond length using three different basis set for 
the three-four-body Jastrow term. 







VMC 


LRDMC 


JHF 
JHF 
JHF 


(34BJ 2s2p) 
(34BJ 3s2p) 
(34BJ 3s3p) 


-11.0239(3) 
-11.0245(2) 
-11.0250(2) 


-11.0539(3) 
-11.0543(2) 
-11.0550(3) 


JAGPn = 
JAGPn = 
JAGPn = 


8 (34BJ 2s2p) 
8 (34BJ 3s2p) 
8 (34BJ 3s3p) 


-11.0535(1) 
-11.0536(2) 
-11.0544(1) 


-11.0748(2) 
-11.0745(1) 
-11.0745(2) 



mental internuclear distance. Total energies are reported 
in Table ED 



b. C 2 convergence in the Jastrow basis set 

We further checked the effects of a larger three and 
four-body Jastrow factor (34BJ) in the case of C 2 . 
We performed simulations at the experimental bond 
length with the JHF wave function and the JAGP with 
n = 8 molecular orbitals (whose results agrees with the 
JAGPn* as shown in Appendix IB 1[) with the 2s2p Jas- 
trow used for all the other cases and for two larger ba- 
sis sets (namely 3s2p, and 3s3p). Results are shown in 
Table IVIIl As expected, effects on the molecular total 
energy of a larger three- and four-body Jastrow are neg- 
ligible in the LRDMC at least for the RVB wave func- 
tion. For all the other cases, the largest Jastrow basis 
provides an energy gain of at most 1 mH with respect to 
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the smaller basis. 
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